Shrinking Alpine chamois: Higher spring temperatures over the last 27 years in Southern Switzerland are linked to a 3 kg reduction in body mass of yearling chamois
Supplementary materials and codes for the manuscript
^1 Department of Biology, University of Ottawa, Canada
^2 Swiss Ornithological Institute, Seerose 1, 6204 Sempach, Switzerland
^3 School of Biological Sciences, University of Aberdeen, United Kingdom
^4 Department of Biology, University of Fribourg, Chemin du Musée 10, Fribourg, Switzerland
^5 Studio alpino Tettamanti, La Campagna d Zora 15, 6678 Lodano, Switzerland
^6 Ufficio della Caccia e della Pesca del Cantone Ticino, Bellinzona, Switzerland (old address)
__*Corresponding author__:
1.1 ORCIDs
GM: 0000-0003-4429-7726
LFB: 0000-0001-9552-8032
PB: 0000-0002-6759-4371
2 Abstract and keywords
Abstract: Although climate change is considered to be partly responsible for the size change observed in numerous species, the relevance of this hypothesis for the ungulates remains debated. We used body mass measurements of 5635 yearlings (i.e. 1.5 years old) Alpine chamois (Rupicapra rupicapra) harvested in September in the Swiss Alps (Ticino canton) from 1992 to 2018. In our study area, during this period, yearlings shrank by ca. 3 kg while temperatures between May and July rose by 1.7°C. We identified that warmer temperatures during birth and the early suckling period (May 9 to July 2 in the year of birth) had the strongest impact on yearling mass. Further analyses of year-detrended mass and temperature data indicate that this result was not simply due to changes in both variables over years, but that increases in temperature during this particularly sensitive time window for development and growth are responsible for the decrease in body mass of yearling chamois. Altogether, our results suggest that rising temperatures in the Alpine regions could significantly affect the ecology and evolution of this wild ungulate.
Keywords: climate change, climwin, ungulates, life stages, temperature, elevation
Journal: Royal Society Open Science
3 Libraries and datasets
3.1 Libraries
knitr::opts_chunk$set(
fig.path = "figures/",
dev = c("png", "tiff", "postscript", "pdf"), # for papers ("png", "tiff")
dpi = 300
)
# load the packages
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(snakecase)
library(climwin)
## Loading required package: ggplot2
## Use suppressPackageStartupMessages() to eliminate package startup messages
## Loading required package: gridExtra
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
## Loading required package: Matrix
## To learn how to use climwin see our vignette.
## See help documentation and release notes for details on changes.
library(tidyr)
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
library(ggplot2)
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(lme4)
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(stringr)
library(MuMIn)
3.2 Session information
R session information is printed here for repeatability.
The data analysed in this study are the records of the Ticino hunting bags from 1992 to 2018. In Ticino, hunting starts at the beginning of September and the harvest plan is mostly completed within three weeks.
Data were collected from the Alps in Ticino, the southernmost canton of Switzerland, over an area of 2700 km2 with an elevation varying from 250 to 2700 m asl. The climate in the mountain range is Alpine, with temperatures varying from mean temperatures of -12°C in winter to mean temperatures of 15.5°C in summer. The hottest and the sunniest month of the year is July with an average maximum temperature of 25°C, measured in the biggest city in the canton Lugano (World Weather & Climate Information, 2021).
Overall, 34 017 animals were legally shot during the hunting period ranging from an age of 0.5 to 22.5 years old. All animals were sexed, aged and weighted (eviscerated). Both males and females have horns all year-round, even though female ones tend to be shorter. For the estimation of the age of the shot chamois, measurement of the teeth and the growth rings of their horns were used (Schroder and Elsner-Schack 1985).
# load the data
ch_biom <- read.csv("data_chamois_yearlings.csv", stringsAsFactors = TRUE, na = c("", "NA"))
clim <- read.csv("data_swiss_weather.csv", stringsAsFactors = TRUE, na = c("", "NA", "-"))
colnames(ch_biom) <- snakecase::to_snake_case(colnames(ch_biom))
# fixing some variables
ch_biom$date_ymd <- as.Date(paste(ch_biom$year, ch_biom$month, ch_biom$day), "%Y %m %d")
clim$date_ymd <- as.Date(clim$date, "%d/%m/%y")
ch_biom$year_f <- as.factor(ch_biom$year)
3.3.1 Subset
Due to the nature of the dataset, only information on individuals shot in September was available, so for the purpose of this study, only a 1.5-year-old animals were considered (7127 individuals, 3257 females and 3870 males). As chamois are usually weaned at 3 to 6 months of age (Scornavacca et al. 2018), a 1.5-year-old individual has been feeding on their own for nearly a year, is fully grown but still very vulnerable to external abiotic and biotic threats due to the decrease in maternal care and increase in active grazing behaviour.
Daily mean ambient temperature (°C) from 1990 until 2018 (all the years needed for the analysis) was obtained from a Swiss meteorological station in the city of Lugano (273 m asl), within the harvesting area.
As this weather station is at a lower elevation compared to the harvesting area of the Chamois, we tested here for correlations with 2 higher elevation stations, both located close to the town of Acquarossa.
The first one is located in Comprovasco (Coordinates: 714984/146451, Elevation: 575m a.s.l).
As both weather station present high correlation values with the station of Lugano, we decided to use this last weather station in the models as it includes all the years necessary for the analyses
##
## Pearson's product-moment correlation
##
## data: clim$temp_mean_lugano and clim$temp_min_lugano
## t = 632.13, df = 18626, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9768293 0.9781087
## sample estimates:
## cor
## 0.9774779
plot(clim$temp_mean_lugano, clim$temp_min_lugano)
5 Supplementary Material 2
As the use of arbitrary climate periods do not always explain the biological response in the best way possible (van de Pol et al. 2016), we investigated the variation weight of yearling individuals in relation to the variation of mean ambient temperature using the R package climwin, and the function slidingwin which detects the exact time window when a biological variable is most strongly affected by climate (Bailey and van de Pol 2016).
The overall approach for the climwin analysis is to compare the support by the data for competing hypotheses and to formalize them into regression models (van de Pol et al., 2016).
Competing models are based upon a baseline model (called also null model, a model without weather effects) and ranked using the deltaAICc, or the difference in terms of the Akaike Information Criterion values calculated for a small sample size between the candidate model and baseline model.
Climwin presents the models using the deltaAICc value relative to the baseline model (AICc of the candidate model - AICc of the baseline model). Therefore, a model that is more supported than the baseline model will have a negative deltaAICc value. On the same hand the model with the best support from the data, usually with lowest AICc, will be shown as the model with lowest deltaAICc in the climwin output.
The baseline model was a linear model with the body mass of the yearling chamois in relation to sex and elevation. The function slidingwin creates a candidate set of competing models testing windows of different lengths for the weather variable of interest, in this study the mean daily ambient temperature (°C).
Non-linear effects of temperature on body weight were taken into account by checking for both linear and quadratic trends. This is mentioned in the climwin output as func = lin (only linear term) func = quad (linear and quadratic terms).
As most of the chamois was shot during a two-week period at the end of September we chose an absolute time window for the analyses instead of an individual specific time window. As reference day we chose the last date of the shooting period (September 24th) and we looked for windows between September 24th and 662 days before (December 1st of 2 years before) to include the three critical periods of a young chamois life: gestation, lactation and yearling.
5.1 Base model
According to (van de Pol et al. 2016), we built a base model that includes variables that can affect the body size, i.e. elevation and sex.
ch_basemod <- lm(
weight ~
sex + elevation_sc,
data = ch_biom15
)
summary(ch_basemod)
##
## Call:
## lm(formula = weight ~ sex + elevation_sc, data = ch_biom15)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.2112 -1.8462 0.0268 1.7888 13.1538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.90777 0.05335 260.668 < 2e-16 ***
## sexm 0.56032 0.07145 7.842 5.28e-15 ***
## elevation_sc 0.47942 0.03549 13.509 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.662 on 5632 degrees of freedom
## Multiple R-squared: 0.04296, Adjusted R-squared: 0.04262
## F-statistic: 126.4 on 2 and 5632 DF, p-value: < 2.2e-16
5.2 Climwin analysis
5.2.1 Finding the best window
Using the function slidingwin allows to search for the best climatic window
The linear+quadratic term better explains the variation in the data (deltaAICc has the lowest value), sorted by deltaAICc such that the best supported model is on top.
To investigate any other tested hypothesis we can simply replace the number in the double square brackets with the corresponding list number.
ch_mass_sw$combos %>% arrange(DeltaAICc)
AICc of the Best model with the linear+quadratic term
MuMIn::AICc(ch_mass_sw[[2]]$BestModel)
## [1] 26704.2
AICc of the Best model with the linear term
MuMIn::AICc(ch_mass_sw[[1]]$BestModel)
## [1] 26767.51
AICc of the baseline model (no climatic factor), used by the function slidingwin as a reference to obtain the deltaAICc values plotted above:
MuMIn::AICc(ch_basemod)
## [1] 27029.52
Difference in terms of AICc between the Best model and the baseline model
DeltaAICc as obtained using the function slidingwin in the climwin package
ch_mass_sw[[2]]$Dataset$deltaAICc[1]
## [1] -325.3275
They are the same!
5.2.4 The 30 best quadratic models
The 30 best windows for the linear+quadratic models sorted by deltaAICc. All models with the lowest AICc (delta AICc between -325.3275 and -320.4684) present very comparable windows: - WindowOpen and WindowClose similar (+- 3 days) to the one of the top model.
head(ch_mass_sw[[2]]$Dataset, 30)
5.2.4.1 Windows plot
It's possible to extract the time windows of all the best supported models (i.e. multi-model inference). This panel shows the opening and closing points of the time windows that were best supported by the data, here those models that made up 95% model confidence set.
plotwin(ch_mass_sw[[2]]$Dataset)
5.2.4.2 Delta plot
The variation in deltaAICc between time windows can be better investigated using the following plot:
Warmer areas shows values with the lowest deltaAICc (i.e. "best models"). As explained by van de Pol et al., 2016, these deltaAICc landscapes of the different time windows shows multiple peaks (red areas) instead of a clear single peak. This can indicate the presence of multiple (e.g. possibly both long- and short-lag) weather signals within the same weather variable, but it can also occur due to collinearity or chance.
The evidence for multiple signals can be therefore investigated by adding the best supported of the weather windows to the baseline model, and re-fitting all the different time windows again: this tests whether there is still strong model support for the second best (e.g. short-lag) weather window once the other best supported (e.g. long-lag) weather window has been accounted for in the baseline model (here in the Step 2).
5.2.4.3 Beta plot
This panel shows the model support (deltaAICc) for all fitted time windows tried, shown for each combination of Window open (y-axis) and Window close (x-axis). Models with the lowest deltaAICc (red) are the best supported (colours show the deltaAICc levels compared to the null model, see legend). Strongly supported windows will often be grouped together.
plotbetas(ch_mass_sw[[2]]$Dataset, arrow = TRUE)
5.2.4.4 Autocollinearity
Correlation between the mean temperature during the best supported time window and the mean temperature over all other time windows.
load(file = "climwin_autocall_tmean.rda")
plotcor(autocoll, type = "A", arrow = TRUE)
5.3 Main results
5.3.1 The best window
Dates of the best window (as if compared to year of harvest 2018)
as.Date("2018/09/24", format = "%Y/%m/%d") - ch_mass_sw$combos$WindowOpen[[2]]
## [1] "2017-05-09"
as.Date("2018/09/24", format = "%Y/%m/%d") - ch_mass_sw$combos$WindowClose[[2]]
## [1] "2017-07-02"
5.3.2 The model
I can add the new temperature variable for the extracted time window to the original dataset:
# The best supported climate variable can be attached
# to the original dataset for further analyses
ch_biom15$temp_503_449 <- ch_mass_sw[[2]]$BestModelData$climate
Relationship between body mass of harvested 1.5-year-old Alpine chamois and (a) the average temperature between May 9 and July 2 of the birth year (suckling period), and (b) elevation (m a.s.l.). Each dot is one observation (darker dots represent a higher number of observations); fitted lines in (a) and (b) are shown with 95 % confidence intervals (shaded areas).
Note that the quadratic model is heuristic and does not imply that the relationship is parabolic over the whole range of temperatures.
5.3.3 Last step: Randwin
Using randwin to randomize the identity of the chamois we are able to check if the window that was found before is actually important, or the relationship was just random.
# Performing randamization to identify
# likelyhood of of dignals occuring by chance
ch_mass_rand100 <- randwin(
repeats = 100,
baseline = ch_basemod,
xvar = list(Temp = clim$temp_mean_lugano),
type = "absolute",
refday = c(24, 9),
range = c(662, 0),
stat = "mean",
cdate = clim$date_ymd,
bdate = ch_biom15$date_ymd,
func = c("lin", "quad"),
cmissing = FALSE,
cinterval = "day",
window = "sliding"
)
save(ch_mass_rand100, file = "climwin_mass_randomization.rda")
## Warning in pvalue(datasetrand = ch_mass_rand100[[2]], dataset = ch_mass_sw[[2]]$Dataset, : Pc will
## be overly conservative when sample size is greater than 47
## [1] 3.144653e-06
The randomization process shows that the window is actually important.
Here we ran year-detrended analyses to demonstrate that year is not confounding the relationship between body mass and temperature. We extracted the residuals of linear regressions between mass and year and between temperature and year, and then ran a linear model with the residuals of body mass in relation to the residuals of temperature.
Annual trend of (a) average temperature between May 9 and July 2 and (b) body mass of harvested 1.5-year-old Alpine chamois between 1992 and 2018, and (c) year-detrended relationship between body mass and temperature. Detrended values in (c) are residuals from linear models in (a) and (b). Each dot is one observation (darker dots representing a higher number of observations in (b)); fitted lines are shown with 95% confidence intervals (shaded areas).
6 Supplementary Material 3
Analyses with the minimum and maximum temperature, same base model as in Supplementary Material 2
6.1 Climwin analysis
6.1.1 Finding the best windows
Using the function slidingwin allows to search for the best climatic window
When considering mean, minimum or maximum temperature, the linear+quadratic term better explains the variation in the data (deltaAICc has the lowest value), sorted by deltaAICc such that the best supported model is on top.
To investigate any other tested hypothesis we can simply replace the number in the double square brackets with the corresponding list number.
ch_mass_sw$combos %>% arrange(DeltaAICc)
6.1.4 Maximum temperature
Dates of this window (as if compared to year of harvest 2018)
as.Date("2018/09/24", format = "%Y/%m/%d") - ch_mass_sw$combos$WindowOpen[[5]]
## [1] "2017-04-20"
as.Date("2018/09/24", format = "%Y/%m/%d") - ch_mass_sw$combos$WindowClose[[5]]
## [1] "2017-07-11"
The maximum temperature has a wider window (earlier Open date and later Close date) compared to the mean temperature, but the window overlaps.
windows plot
It's possible to extract the time windows of all the best supported models (i.e. multi-model inference). This panel shows the opening and closing points of the time windows that were best supported by the data, here those models that made up 95% model confidence set.
Interpretation: Warmer areas shows values with the lowest deltaAICc (i.e. "best models"). As explained by van de Pol et al., 2016, these deltaAICc landscapes of the different time windows shows multiple peaks (red areas) instead of a clear single peak. This can indicate the presence of multiple (e.g. possibly both long- and short-lag) weather signals within the same weather variable, but it can also occur due to collinearity or chance.
The evidence for multiple signals can be therefore investigated by adding the best supported of the weather windows to the baseline model, and re-fitting all the different time windows again: this tests whether there is still strong model support for the second best (e.g. short-lag) weather window once the other best supported (e.g. long-lag) weather window has been accounted for in the baseline model (here in the Step 2).
Dates of this window (as if compared to year of harvest 2018)
as.Date("2018/09/24", format = "%Y/%m/%d") - ch_mass_sw$combos$WindowOpen[[6]]
## [1] "2017-05-19"
as.Date("2018/09/24", format = "%Y/%m/%d") - ch_mass_sw$combos$WindowClose[[6]]
## [1] "2017-06-30"
The maximum temperature has a narrower window (laterer Open date and earlier Close date) compared to the mean temperature, but the window overlaps.
windows plot
It's possible to extract the time windows of all the best supported models (i.e. multi-model inference). This panel shows the opening and closing points of the time windows that were best supported by the data, here those models that made up 95% model confidence set.
Interpretation: Warmer areas shows values with the lowest deltaAICc (i.e. "best models"). As explained by van de Pol et al., 2016, these deltaAICc landscapes of the different time windows shows multiple peaks (red areas) instead of a clear single peak. This can indicate the presence of multiple (e.g. possibly both long- and short-lag) weather signals within the same weather variable, but it can also occur due to collinearity or chance.
The evidence for multiple signals can be therefore investigated by adding the best supported of the weather windows to the baseline model, and re-fitting all the different time windows again: this tests whether there is still strong model support for the second best (e.g. short-lag) weather window once the other best supported (e.g. long-lag) weather window has been accounted for in the baseline model (here in the Step 2).
We thank the managers of the hunting and fishing cantonal office of Ticino, Switzerland, and the Swiss federal office of meteorology and climatology (MeteoSchweiz) for collecting the data and making them available to us.
8 Funding
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101025938 to GM.
9 Data accessibility
All data and code used for statistical analysis and plots are provided via the Open Science Framework at https://osf.io/p9c4m/ and were shared with editor and reviewers at first submission.
10 Authors' contributions
G.M. and P.B. conceived the study. F.T. compiled the data, and L.F.B and N.I curated the data. G.M. and K.G.G. performed the statistical analyses with the help of P.B. G.M. and K.G.G. drafted the manuscript, and all authors provided inputs at all stages. All authors approved the final version of this manuscript, and all authors agree to be held accountable for the work performed therein.
